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Abstract 

We consider the random-bond it J Ising model on a square lattice as a function of the temperature 
T and of the disorder parameter p {p = 1 corresponds to the pure Ising model). We investigate 
the critical behavior along the paramagnetic-ferromagnetic transition line at low temperatures, 
below the temperature of the multicritical Nishimori point at T* = 0.9527(1), p* = 0.89083(3). We 
present finite-size scaling analyses of Monte Carlo results at two temperature values, T ~ 0.645 and 
T = 0.5. The results show that the paramagnetic- ferromagnetic transition line is reentrant for T < 
T*, that the transitions are continuous and controlled by a strong-disorder fixed point with critical 
exponents = 1.50(4), rj = 0.128(8), and (3 = 0.095(5). This fixed point is definitely different 
from the Ising fixed point controlling the paramagnetic- ferromagnetic transitions for T > T*. Our 
results for the critical exponents are consistent with the hyperscaling relation 2(3 /v — r] = d — 2 = 0. 

PACS numbers: 75.10.Nr, 64.60.Fr, 75.40. Cx, 75.40.Mg 
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I. INTRODUCTION 



The ±J Ising model represents an interesting theoretical laboratory, in which one can 
study the effects of quenched disorder and frustration on the critical behavior of spin systems. 
While originally introduced to describe magnetic systems with disordered couplings [1], it 
has been shown recently to be also relevant for quantum computations [2, 3]. It is defined 
by the lattice Hamiltonian [1] 

7i = — ^ Ja;yCra;CTy, (l) 

where — ±1, the sum is over all pairs of lattice nearest-neighbor sites, and the exchange 
interactions J^y are uncorrelated quenched random variables, taking values ± J with proba- 
bility distribution 

P{J.y) = pS{J.y - J) + (1 - p)S{J,y + J). (2) 

In the following we set J = 1 without loss of generality. For p = 1 we recover the standard 
Ising model, while for p = 1/2 we obtain the bimodal Ising spin-glass model. 

The T-p phase diagram of the two-dimensional (2D) square-lattice ±J Ising model has 
been extensively investigated [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. 
The resulting phase diagram, which is sketched in Fig. 1, presents two phases at finite 
temperature: a paramagnetic and a ferromagnetic phase. They are separated by a transition 
line, which starts at the pure Ising transition point at p = 1 and Tjs ~ 2.269 and ends at 
the T = transition at Pq ^ 0.897. The point where this transition line meets the so-called 
Nishimori (N) line [13], at T* = 0.9527(1) and p* = 0.89083(3) (we derive these estimates 
in the present paper), is a multicritical point (MNP) [36]. 

The MNP divides the paramagnetic-ferromagnetic (PF) transition line in two parts. The 
PF transition line from the Ising point at p = 1 to the MNP is controlled by the Ising 
fixed point: disorder gives only rise to logarithmic corrections to the standard Ising critical 
behavior [11]. On the other hand, the presence of the MNP on the transition line suggests 
that the PF transitions for T < T* belong to a different strong-disorder universality class. 
This is confirmed by the renormalization-group (RG) calculations of Refs. 41 and 30, using 
domain-wall and Migdal-Kadanoff RG transformations respectively, which found that the 
RG flow along the critical line for T < T* was attracted by a different fixed point. 
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FIG. 1: (Color online) Phase diagram of the square-lattice it J Ising model for 1 ~ p < 1/2. 

In this paper we investigate the critical behavior along the low-temperature transition 
line from the MNP to the T = axis. We perform Monte Carlo (MC) simulations at two 
temperature values below the MNP, i.e., ai (3 = 1/T = 2 and (3 = 1.55. As we shall see, our 
finite-size scaling (FSS) analyses show that the PF transition line for T < T* is reentrant 
and that the transitions are continuous. Moreover, the estimates of the critical exponents 
and of several RG invariant quantities for these two values of T are consistent, supporting 
the hypothesis that the PF transition line below the MNP belongs to a unique universality 
class. The values of the critical exponents, u = 1.50(4), r] = 0.128(8), and P = 0.095(5) are 
clearly different from the Ising values z/ = 1, 77 = 1/4, /? = 1/8. Therefore, these results 
show the existence of a strong-disorder fixed point associated with a PF transition. Note 
that this strong-disorder fixed point does not violate hyperscaling. Indeed, our results are 
consistent with the hyperscaling relation 2 + 2f3/v — iri = d = 2 (our estimates of the critical 
exponents 7] and (3 give 2 + 2/5/z/ — 77 = 2.00(1)). The transitions for T <T* are no longer in 
the basin of attraction of the Ising fixed point, which is the relevant one for small disorder 
and determines the critical behavior along the transition line for T > T*. 

The paper is organized as follows. In Sec. II we review the main features of the T-p phase 
diagram of the square-lattice ±J Ising model. The MC results and their FSS analyses are 
presented in Sec. III. In Sec. IV we draw our conclusions. Some technical details on the 
simulations are presented in App. A, while the quantities we compute are defined in App. B. 
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In App. C we present a reanalysis of the critical behavior at the MNP, using the additional 
data we have collected in this work. Moreover, we also present analyses which take into 
account the analytic corrections, which had been neglected in our previous work [15]. This 
allows us to obtain improved estimates of the critical parameters at the MNP. 



II. THE PHASE DIAGRAM OF THE SQUARE-LATTICE ±J ISING MODEL 

The phase diagram of the square-lattice ±J Ising model is sketched in Fig. 1. It is 
symmetric for p — > 1 — p and thus we only report it for 1 — p < 1/2. For sufficiently 
small values of the probability of antiferromagnetic bonds pa = 1 — p, the model presents 
a paramagnetic phase and a ferromagnetic phase, separated by a transition line. The PF 
transition line starts at the Ising point Xjs — {T — Tis,p— 1), where Tis = 2/ ln(l + ^/2) — 
2.26919... is the critical temperature of the 2D Ising model, and extends up to a T = 
transition at [42, 43] Xq = {T = 0,p = po ^ 0.897). 

The slope of the transition line at p = 1 is known exactly [45] , so that for small 1 — p we 
have 



1 ^(l-p) + 

ln(H-V2)^ 



(3) 



In the T-p phase diagram an important role is played by the Nishimori (N) line [4, 13] 
defined by the equation (p > 1/2) 

T = r.(p), r.(p) = j-— A_. (4) 

Along the N-line several rigorous results can be proved [4, 13, 46]. The energy density is 
given by 

EN{p) = ^[{n)T^ip)]^2-ip, (5) 
and the spin-spin and the overlap correlation functions are equal 

[(o-qo-x)] = [{(JqCTx)^]- (6) 

Here the angular and square brackets refer respectively to the thermal average and to the 
quenched average over the bond couplings {Jxy}- As argued in Refs. [36, 39] and verified 
numerically [15, 18, 19, 20, 22, 33], the critical point Xmnp = {T* ^ 0.953, p* ^ 0.891) along 
the N line is a multicritical point (MNP). 
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Along the transition line from the Ising point Xis to the MNP, the critical behavior is 
analogous to that observed in 2D randomly dilute Ising (RDI) models [11]. It is controlled 
by the pure Ising fixed point and disorder is marginally irrelevant, giving rise to a universal 
pattern of logarithmic corrections, see, e.g., Refs. [11, 47, 48, 49] and references therein. 

The location of the MNP and the corresponding critical exponents can be obtained by 
FSS analyses of MC data along the N line. The new analysis reported in App. C gives 

T* = 0.9527(1), p* = 0.89083(3). (7) 

In the absence of external fields, the MNP is characterized by two relevant RG operators 
with RG dimensions yi — 0.66(1) and y2 = 0.250(2). Moreover, the magnetic exponent r] is 
given hj 1] = 0.177(2). Other estimates of T*, p*, and of the critical exponents can be found 
in Refs. [14, 15, 16, 18, 22]. 

As a consequence of the inequality [13] 

\[{<^x<^y)T]p\ < [\{<^x<^y)TMip)\]p (8) 

(the subscripts indicate the vahics of T and p at which the thermal and disorder average are 
performed), ferromagnctism can only exist in the region p > p*. Thus, the PF boundary 
lies in the region p > p* and, at the MNP, the transition line is tangent to the line p — p*, 
hence parallel to the T axis. As a further consequence, at T = the ferromagnetic phase 
ends at p = po with po > p*. In Refs. [4, 34, 35, 38] it was argued that the PF transition 
fine from the MNP to Xq — (0, po) is only related to the frustration distribution; hence, it 
should not depend on temperature and should coincide with the line p = p*, so that po = p*. 
Numerical estimates of Pq have shown that this argument is not exact. Indeed, numerical 
analyses [12, 18, 22, 42, 43, 44] give po ~ 0.897;^ this suggests that the transition hne below 
the MNP is reentrant, i.e. Pc > p* for any T < T*. The difference is however quite small, 
Po - P* ~ 0.006. 

Our FSS analyses confirm that the PF transition line is reentrant for T < T*. Indeed, 
we find pe = 0.8915(2) at T = 1/1.55 ^ 0.645 and Pe = 0.8925(1) at T = 0.5. The PF 

transitions are of second order and show the same critical behavior with critical exponents 
u = 1.50(4), r] = 0.128(8), and j3 = 0.095(5), which are consistent with hyperscaling. These 

^ The most precise estimates are apparently [42] po = 0.897(1) and [43] po = 0.8969(1). 
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results confirm the existence of a strong-disorder fixed point, different from tlie Ising fixed 
point which controls the PF transitions above the MNP, i.e. for T* <T < Tig. 

At variance with the three-dimensional case, there is no evidence of a finite-temperature 
glassy phase. Glassy behavior is only expected for T = and p < po- The critical behavior 
for T — > has been much investigated for p = 1/2 [6, 7, 8, 9, 10]. In particular, simulations 
found that the correlation length increases as T~'^ with u ^ 3.5. A natural hypothesis is 
that a T = glassy transition occurs for any p < po, with critical behavior in the same 
universality class as that of the bimodal model with p — 1/2. 

The point Xq — (0,po)) where the low-temperature transition line ends is a multicritical 
point: it is connected to three phases and it is the intersection of two different transition 
lines, the PF line at T > and the glassy hue at T = 0. At T = the critical point Xq 
separates a ferromagnetic phase from a T = glassy phase, while for T > the transition 
line separates a ferromagnetic from a paramagnetic phase. Therefore, on general grounds, 
the critical behavior when varying p at T = differs from that along the PF transition 
line at T > 0, unless the magnetic and glassy critical modes are effectively decoupled at 
the T — multicritical point. The latter scenario is apparently supported by the fact that 
the estimates of magnetic critical exponents at T = 0, see e.g. Refs. 18, 41, 42, 43, are 
quite close and substantially consistent with those found along the transition line at finite 
temperature < T < T* . 

III. MONTE CARLO RESULTS 

We investigate the critical behavior along the PF fine that starts at the MNP T* 0.95 
and ends at T = 0. Since the transition fine below the MNP is expected to be almost 
parallel to the T axis, we study the FSS behavior of several quantities at fixed T as a func- 
tion of p. We consider two values of T, j3 = 1/T = 2 and f3 = 1.55, which are quite far 
from the two endpoints of the line. For each of these two values we perform MC simula- 
tions on square lattices of linear size L with periodic boundary conditions, for several values 
of L: L — 8, 12, 16, 24, 32, 48, 64. In our MC simulations we employ the Metropolis algo- 
rithm, the random-exchange method (often called parallel-tempering or multiple Markov- 
chain method) [50, 51], and multispin coding. Some details are reported in App. A. 
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A. The critical point Pc and exponent v 



Wc first focus on the data at /3 = 2, for wliicli we have most of the statistics. In order 
to estimate Pc and i/, we perform a FSS analysis of the renormahzed couphngs i?^, t/4, t/22, 
which are defined in App. B and are generically denoted by R in the following. MC estimates 
are shown in Fig. 2. We clearly observe a crossing point for 0.8920 <p < 0.8930, indicating 
Pc ^ 0.892-0.893, which is larger than the value at the MNP, i.e. = p* = 0.89083(3). This 
already suggests that the transition line is reentrant. 

To obtain more precise estimates wc perform a careful FSS analysis, following Ref. [52]. 
In the FSS limit any RG invariant quantity obeys the scaling law 

R = fR{uiLy^), (9) 

where /r(0) = R*, yi = and we have neglected scahng corrections. Here Ui is the 
nonlinear scaling field associated with the leading relevant operator, which has RG dimension 
Hi. The scaling field is an analytic function of the system parameters which vanishes along 
the critical line. Thus, for p Pc{P) at fixed /9 we can write 

ui = Ao{(3){p - Pc) + A,{(3){p - pcf + ... (10) 

where the coefficients A(/9) are analytic functions of (3. The terms of order {p — p^y, {p — 
PcY, etc., give rise to corrections of order L'^^^ as L ^ 00. They arc named analytic 
corrections, because they arise from the analytic dependence of the scaling fields on the model 
parameters. See Ref. [52] for a thorough discussion of their origin. In pure ferromagnetic 
systems, in which v < 1 and yi 1, they are usually negligible, and the nonanalytic 
corrections, which behave as L~'^, c<; < 1, play a much more important role. This is not the 
case here, since, as we shall see, at the transition line ?/i ~ a; < 1. 

Since our data are sufficiently close to the critical point, p — pc'is small and thus wc can 
take ui ~ [p — Pc)- Moreover, also the product (p — Pc)Ly^ is small, so that we can expand 
fR{x) in powers of x. Thus, we fit the numerical data to 

R^R^ + J^anip-PcTL"^'', (11) 

71=1 

keeping R*, the coefficients {a„}, Pc, and yi as free parameters. Here we neglect scaling 
corrections. To monitor their role, we repeat the fits several times, each time only including 
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FIG. 2: (Color online) MC estimates of = i/L, C/4, U22 at /3 = 2 vs p. The lines connecting 
the data at given L are drawn to guide the eye. The dashed vertical lines corresponds to the MNP 
location p* = 0.89083(3). The dotted vertical lines indicate our final estimate of Pc = 0.8925(1). 
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48/54 
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TABLE I: Estimates obtained from the analysis of tlic data at /3 = 2. Above we report the results 
of the combined fits of R^, U4, and U22 to Eq. (11) with Wmax = 2. Below we report the results of 
the fits to Eq. (12) with rimax = 2 and kmax = 1- 

data satisfying L > Lmm- For a given Lmin, X^/DOF (DOF is the number of degrees of 
freedom of the fit) changes significantly as we increase rimax from 1 to 2, and only marginally 
as we change this parameter from 2 to 3. This indicates that the range of values of p we are 
considering is too large to allow for a linear approximation of the scaling function fn^x). 
Instead, a quadratic approximation seems to be accurate enough. Thus, the results we 
present below correspond to Wmax = 2. 

In Table I we give the estimates of R* , Pc, and yi from combined fits of R^, U4, and 
U22- All quantities, except yi, show a significant — much larger than the statistical errors — 
variation with Lmin- Moreover, the is very large. Clearly, scaling corrections are not 
neghgible. In order to take them into account, we fit the MC data to 

R^R* + Y^ a,{p - p^fL^y^ + ^^(^ - PcfL^'"\ (12) 

n=l k=0 

taking a; as a free parameter. Results for fcmax = 1 and rimax = 2 are also reported in Table 
I. The is now significantly smaller and x^/DOF 1, indicating that the fitting form (12) 
describes the data at the level of their statistical accuracy. The results are stable and the 
estimates for Lmin > 12 are consistent within errors. These fits also provide an estimate of 
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the correction-to-scaling exponent lo. We find 



a; = 0.6(1). (13) 

The estimates of lo and y\ indicate that a; ~ so that analytic and nonanalytic corrections 
behave analogously. Therefore, we should also consider the analytic corrections. For this 
purpose, we also performed fits to 

^max 

R = K* + Y. <^ + - ^^^)]"(^^ - VcTL^'>\ (14) 

n=\ 

which corresponds to including the quadratic term in the expansion of the nonlinear scaling 
field u\. The parameter c is a new fitting parameter which is independent of the quantity 
one is analyzing. Fits to Eq. (14) are substantially equivalent to those to Eq. (11). For 
instance, the of the combined fit for Lmin = 8 is 9846, which is identical to that reported 
in Table I for the same value of -Lmin- The coefficient c is small and we estimate |c| < 0.3. 
Since our data satisfy |p — pd < 0.0030, the analytic term gives a tiny correction and does 
not influence the fit results. 

Comparing the results of the different fits we arrive at the final estimates 



2/1 = 0.67(2), V = 1/y, = 1.50(4), (15) 

Pe = 0.8925(1), (16) 

= 1.25(3), (17) 

U* = 1.084(3), (18) 

= 0.074(1). (19) 

The central value corresponds to the result of the fit to Eq. (12) with Lmin = 12; the errors 



are such to include the results of the fits to Eq. (11) and Lmin = 32, and should take into 
account the systematic error due to further scaling corrections which have been neglected 
in our analyses. 

We repeat the same type of analysis at P — 1.55. We report in Table II the results of the 
fits to Eqs. (11) and (12). In the latter case the data do not allow us to perform fits in which 
u; is a free parameter. Thus, wc only report results of fits in which uj is fixed to 0.4, 0.6, and 
0.8, consistently with the estimate a; ~ 0.6 presented above. Fits without scaling corrections 
are characterized by large values of x^/DOF and by a systematic trend of the results. Fits 
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236/43 1.1543(15) 1.0959(3) 0.0755(2) 0.89131(3) 0.623(17) 


8 


0.4 


77/82 1.213(3) 


1.0878(5) 0.0739(4) 0.89160(3) 0.663(52) 
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1.0908(4) 0.0744(3) 0.89153(2) 0.662(37) 
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1.0904(3) 0.0746(2) 0.89148(2) 0.659(29) 



TABLE II: Estimates obtained from the analysis of the data at /3 = 1.55. Above we report the 
results of the combined fits of R^, U4, and U22 to Eq. (11) with rimax = 2. Below we report the 
results of fits to Eq. (12) with nmax = 2, fcmax = 1, and lo fixed to 0.4, 0.6, 0.8. 

with scaling corrections are significantly better. The estimates of yi and f/22 perfect 
agreement with those obtained at = 2. Those of f/4 and are substantially consistent: 
the difference between the estimates (17), (18) and the results of the fit with lu — 0.6 — this 
is the fit which, in principle, should be more reliable — is of the order of two error bars and 
can thus be explained by the presence of residual scaling corrections which are not taken 
into account in our error estimate. Therefore, our analyses of the renormalized couplings 
are consistent with a critical transition line whose nature is T independent: for T <T*, the 
PF transition belongs to a unique universality class. 

The estimate (15) is different from the Ising value z/ = 1. Therefore, the PF fixed 
point associated with the transitions along the line T < T* is a new one, clearly distinct 
from the Ising one, which controls the critical behavior for weak disorder. Analogously, 
our estimates of the critical value of the renormalized couplings differ from the Ising values 
[11, 53] i?| = 0.9050488292(4), Ul = 1.167923(5), f/^s = 0, and from those at the MNP (see 
App. C), which are = 0.997(1), Ul = 1.1264(4), and U;2 = 0.0817(3). 

Our analyses also give an estimate of Pc ior j3 = 1.55: 

Pc = 0.8915(2). (20) 

Therefore, for both values of (3 we find pc > p* = 0.89083(3). Thus, the PF transition hue 
is reentrant, contradicting the conjecture of Refs. [4, 34, 35, 38]. 
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B. The exponent rj 



We determine the critical exponent f] from the critical behavior of the susceptibility X- As 
discussed in Ref. [52] in the context of the three-dimensional paramagnetic-glassy transition, 
close to the critical point the susceptibility x behaves as 

X-^ulL'-^fAiP-Pc)L'^ (21) 

where Uh is a function of p related to the magnetic nonlinear scaling field. Note that we 
have approximated ui with p — pc, because, as already discussed, the analytic dependence 
of the scaling field ui is negligible for our data. 

Since we are very close to the critical point, we can expand all quantities in powers of 
(p ~ Pc) ■ For this reason we perform fits to 

In X = (2 - 77) In L + ^ a„(p - pJ^L^^^ + ^mip - pT ■ (22) 

n=0 m=l 

As before, we first analyze the data at /3 = 2. To understand the role of the analytic 
corrections, we first perform fits of the data in which we fixpc = 0.8925 and yi = 0.67, which 
are the estimates obtained above. If we do not include the analytic correction (we set 6m = 
for any m) and we use Timax = 2, we obtain x^/DOF = 633/55, 297/43 from the analysis 
of the estimates of x corresponding to lattices such that L > Lmm = 8, 12, respectively. If 
instead we include the analytic corrections taking mmax = 1, we obtain x^/DOF = 42/54, 
30/42. The improvement is clearly significant, indicating that the analytic corrections cannot 
be neglected. 

In Table III we report the results of the fits corresponding to rimax = 2 and mmax = 1- 
In all cases we fix tji to 0.67(2), as indicated by Eq. (15). This is not crucial, since the 
estimates of the exponent 7] are quite insensitive to this parameter. The results show instead 
a significant dependence on Pc and thus, we present fits in which p^ is fixed to the value 
(16) and fits in which Pc is a free parameter. The estimates of the two fits are substantially 
consistent and show a tiny dependence on Lmin- Also the estimates of Pc are consistent with 
the value (16). 

We also considered nonanalytic scaling corrections, performing a fit of the form 

'^max ^max 

In X = (2 - 7^) In L + J] a^ip - p.yL^y^ + bm{p - Pc)™ + cL'^- (23) 

n=0 m=l 

12 



Pc = 0.8925(1) Pc free parameter 
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TABLE III: Estimates of rj from fits to Eq. (22) with nmax = 2 and mmax = 1- We fix yi = 0.67(2) 
in both fits. The reported error takes into account the error bar on yi and on Pc (for the fit in 
which this quantity is fixed). Analyses of the data at /3 = 2. 

We fix Pc — 0.8925, yi — 0.67, u; — 0.6, ^max = 2, mmax = 1, and obtain rj — 0.1234(13), 
c — 0.000(3) for Lmin — 8: there is no evidence of nonanalytic scaling corrections. 

To avoid the use of pc, note that Eq. (9) can be inverted to give UiL^^ ~ (p — Pc)Ly^ as 
a function of R. Thus, Eq. (21) can also be rewritten as 

X = ulL'-^g^{R)[l + 0{L--)l (24) 

where it! is a renormalized coupling. A polynomial approximation for Uh{p) and QxiR) gives 
the fitting form 

Inx = (2 - 77) InL + J] a„i?" + ^ bmP"" ■ (25) 

n=0 m=l 

Fits to this form have a quite large x^, which is not unexpected since we already found 
that the renormalized couplings show significant scaling corrections. Moreover, the results 
depend significantly on the minimum lattice size Lmin of the data included in the fit. Scaling 
corrections must therefore be included. We thus consider 

'^max ''^max ^max 

In X = (2 - r;) In L + ^ a„i?" + bmP"" + L"" YI (^6) 

n=0 m=l fe=0 

The results of these fits are reported in Table IV. The is good; moreover, the results 
do not depend on which quantity is used in the fit, are stable with Lmm, and are consistent 
with those reported in Table III. 

Analogous analyses can be performed at (3 = 1.55. Also in this case the analytic correc- 
tions cannot be neglected and thus we only consider fits with mmax = 1- The results of the 
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xVDOF 77 


xVdof 7/ 


8 


72/52 0.1247(5) 


47/52 0.1245(14) 


12 


57/40 0.1250(8) 


40/40 0.1253(22) 


16 


51/28 0.1255(9) 


36/28 0.1255(48) 



TABLE IV: Estimates of rj from fits to Eq. (26) with n-max = 2, mmax = Ij ^max = 0, w free 
parameter. On the left we use R = U4, on the right we use R = R^. Analyses of the data at /3 = 2. 



Pc = 0.8915(2) pc free parameter 



'min 


xVdof 


V 


xVdof 


V Pc 


8 


25/28 


0.1336(26) 


23/27 


0.1342(7) 0.89145(3) 


12 


23/23 


0.1335(30) 


22/22 


0.1341(9) 0.89145(4) 


16 


21/18 


0.1335(33) 


20/17 


0.1341(12) 0.89146(6) 


24 


16/13 


0.1330(38) 


16/12 


0.1315(23) 0.89158(12) 


32 


10/8 


0.1330(42) 


10/7 


0.1361(46) 0.89133(21) 



TABLE V: Estimates of t] from fits to Eq. (22) with Ji-max — 2 and ttt-hieix — 

1. We fix yi = 0.67(2) 

in both fits. The reported error takes into account the error bar on yi and on pc (for the fit in 
which this quantity is fixed). Analyses of the data at /3 = 1.55. 

fits to Eq. (22) are reported in Table V. The dependence of the results on Lmin is tiny. More- 
over, the estimates of Pc obtained in the analyses in which this quantity is a free parameter 
are perfectly consistent with the estimate (20). Similar, though less stable, results are ob- 
tained by fitting the data to Eq. (26). We fix a; = 0.6(1) as in the case of the analyses of the 
rcnormalized couplings. For Lmin = 8 we obtain 77 = 0.1304(5) and f] = 0.1340(5) by using 
U4 and R^, respectively; for Lmm = 12 we obtain instead 1] = 0.1317(8) and rj = 0.1329(8) . 

Collecting all results, from the analyses of the data at = 2 we would estimate 77 = 
0.125(3). The analyses at /3 = 1.55 give a slightly different value, rj = 0.132(4). The 
difference is tiny — less than two combined error bars — but indicates that there are corrections 
which are not fully taken into account by our analyses. As final estimate we report the 
average of the two results, 

77 = 0.128(8). (27) 
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The error we quote is quite conservative and essentially includes the estimates of all fits for 
both values of (3. 



C. The exponent /? and a check of hyperscaling 

The exponent f3 can be determined from the critical behavior of the magnetization. The 
RG predicts 

m = UhL-^/-'fm[{p-Pcm, (28) 

where Uh is the same function which appears in Eq. (21) and fm{x) is a universal function. 
Expanding this scaling relation around the critical point we obtain the fitting form 

lnm = -^lnL+^ a„ip - p^^L^y^ + J] 6^(p - p^ . (29) 

n=0 m=l 

As before we fix nmax = 2, mmax = 1, and use the best available estimates of Pc- For (5 — 2 
a fit of the data satisfying L > Lmin — 24 gives P/u — 0.0613(11); for Lmm = 32 we obtain 
instead f3/i^ = 0.0614(12). For f3 = 1.55 and L^in = 32 we obtain = 0.0661(22). As in 
the case of r], we observe a tiny difference between the estimates obtained at the two values of 
the temperature. It probably indicates the presence of additional scaling corrections which 
are not taken into account by our scaling Ansatz. A conservative estimate of the critical 
exponent which is consistent with all results is 

^ = 0.063(3) (3 = 0.095(5). (30) 

We can now check hyperscaling. If it holds, we should have 2(3/i' — ri + 2 = d = 2. We find 

23 

-jy + 2 = 2.00(1). (31) 

Hyperscaling is verified quite precisely. 

Finally, we consider the specific heat. At p = Pc we expect 

a = a + 6L"/^ (32) 

where a is due to the analytic contribution to the free energy. If hyperscaling holds, we 
should have a — 2 — 2u, so that 

^ = ^-2 = 2^1-2 = -0.66(4). (33) 
15 
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0.02 0.04 0.06 0.08 0.1 0.12 0.14 

FIG. 3: The specific heat at /3 = 2 and p = 0.8925 ~ pc versus L^'^l^ . The dashed hne is obtained 
by fitting the MC data to a + bL'"^/^. 

A precise determination of a/u from the data is quite difficult, because a/u < — the 
singular part decreases as L ^ oo. Thus, we have only checked that our data are consistent 
with hyperscaling. In Fig. 3 we show the specific heat for (3 = 2 and L > 24 versus 
The results are consistent, supporting hyperscaling. 




D. The derivative dpc/dP 

As a final test of our results we consider the derivative with respect to /5 of a renormalized 
coupling R. In the FSS limit R behaves as 

R = fniu.Ly^) + u^L-^h,^{u,Ly^) + ... (34) 

where the scaling fields Ui and are functions of the system parameters, hence of (3 and 
p. Moreover, ui vanishes on the critical line. From Eq. (34) we obtain 

H = [fki^^Ly^) + u.L^'' f'nJu.Ly^)] + ^-^L-fnAuiiy^) + ■■■ (35) 

If the critical value pc is /3 independent, duijdp vanishes on the critical line, so that dR/dj3 
behaves as L^'^ for L oo, i.e. the derivative vanishes in the critical large-L limit. This is 
not surprising, since for p = Pc = p* and any (3 we would have R = R* + 0{L^^)^ with R* 
independent of j3. On the other hand, if the transition is reentrant, dR/d(3 diverges as 

We have checked the validity of Eq. (35) by using the data aX. (3 = 2. The fits of the 
renormalized couplings R give us estimates of the expansion of R around Pc- In particular, 
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fits to Eq. (12) give us estimates of the coefficients a„. We have thus fitted dR/dj3 to the 
following expression: 

— = koiy^iai + 2a2(p - Pc)Ly^] + hiy^-^. (36) 

We take yi = 0.67(2), = 0.8925(1), uj = 0.6(1), and ai and as from the fits of R to 
Eq. (12); /cq and ki are free parameters. The estimates of /cq do not vary significantly 
with Lmin. Moreover, results obtained by using dR^/dp and dUi/d(3 are fully consistent. 
Comparing all results we obtain the estimate 

kQ = -0.0020(3) . (37) 

To interpret this result, note that Eqs. (9), (10), and (12) allow us to identify 

ai = A,{(3)f'^{Q). (38) 

Instead, comparing Eq. (36) with Eq. (35) we obtain 

dui 



koai = 



Now, Eq. (10) gives 



d(3 

dui 



/«(0). (39) 

Pc 



^=-MP)%. (40) 



It follows 



dl3 

^ = -A;o = 0.0020(3) . (41) 

Again, this result shows that the transition is reentrant. It is also consistent with the crude 
estimate 

IV. CONCLUSIONS 

In this paper we have studied the nature of the transition line which starts from the MNP 

and ends at T = and which separates the paramagnetic phase from the ferromagnetic 
phase. For this purpose, we have presented FSS analyses of MC data on lattices of linear 
size L up to L = 64 for /3 = l/T = 2 and /3 = 1.55. 
Our main results are the following. 
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(i) The PF transition line below the MNP is reentrant. Indeed, we find Pc = 0.8915(2) at 
T = 1/1.55 ^ 0.645 and Pc = 0.8925(1) at T = 0.5. Therefore, p^ > p* ^ 0.89083(3) 
for any T < T* ^ 0.9527(1), where X* = {T*,p*) is the location of the MNP. 

(ii) The PF transitions are of second order with a standard power-law behavior. 

(iii) The estimated values of the critical exponents and of the large-L limit of the RG 

invariant quantities 1/4^, R^, and U22 at two difi^erent points of the line {f3 = 1.55 and 
/3 = 2) suggest that the PF transitions for < T < T* belong to a unique universality 
class. In particular, the corresponding critical exponents are 

i/ = 1.50(4), 7^ = 0.128(8), /3 = 0.095(5). (43) 

They satisfy the hyperscaling relation 2/3/1/ — rj — d — 2 — 0. Our MC data are 
also consistent with the hyperscaling relation a = 2 — di' = 2 — 21/, which gives 
a = —1.00(8). Using the scaling relation 7 = {2 — 7])v, we derive 7 = 2.81(8). The 
estimates (43) are definitely different from the Ising values u = 1, r] = 1/4, (3 = 1/8. 
We note that they are consistent with the simple rational expressions i/ = 3/2,77 = l/8. 

(iv) The above results show that in two dimensions there are two fixed points which control 
the PF transitions in disordered random-bond Ising systems: besides the standard 
Ising fixed point, which is relevant for small disorder and controls the critical behavior 
along the PF transition fine for T* < T < Tjs, there is also a strong-disorder fixed 
point which controls the critical behavior along the PF transition line for < T < T*. 
The resulting phase diagram is consistent with the results of Refs. [30, 41]. Note that 
frustration and not simply disorder is the relevant property, which gives rise to the 
new fixed point. Indeed, in randomly-dilute Ising systems, in which there is dilution 
but not frustration, there is no evidence of a new strong-disorder fixed point [11]. 

It is interesting to compare our results with those obtained at T = 0. Ref. [41] extrapo- 
lated the RG results to T = (this is correct under the assumption that the limit T — 
is regular) and obtained u — 1.42(8). Ref. [43] obtained u — 1.46(1) from the scaling of 
the failure probability. Ref. [42] found u — 1.55(1) from the analysis of the Binder cumu- 
lant and the magnetization exponent /3 = 0.09(1).^ They also analyzed the domain-wall 

^ Ref. [42] reports /3 = 0.9(1). Alexander Hartmann communicated to us that the correct result is /3 = 
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FIG. 4: The phase diagram of the square-lattice itJ Ising model. The estimates of the critical 
points for X > T* are taken from Refs. [11, 22]. The estimate of the T = transition point is taken 
from Refs. [42, 43]. The dashed lines are interpolations discussed in App. C, while the dotted line 
starting at the Ising point corresponds to the approximation (3). 

energy, obtaining AE = f {{p — pq)L'^^) with p = 0.12(5) and yi = 0.75(5). This gives 
u = l/ui = 1.33(9). The exponent rj associated with the spin-spin correlation has been 
estimated in Ref. [18], obtaining rj ^ 0.13. Note that these estimates are consistent with the 
hyperscaling relation 2(3 /u — ri = d — 2 = 0. Therefore, even if these results refer to a T = 
transition, the magnetic exponents are consistent with hyperscaling. 

The T = results are very close to ours. Note also that, at T = 0, the relation Ul = f/22+l 
holds. ^ This relation is approximately satisfied by our finite-T data, see Eqs. (18) and 
(19); the slight discrepancy might be due to the presence of neglected additional scaling 
corrections. All results are therefore consistent with a single magnetic fixed point that 
controls the magnetic critical behavior both at T > and at T = 0. At the multicritical 
T = point, glassy and magnetic modes are apparently effectively decoupled. 

Finally, we have improved the estimates of the critical parameters at the MNP, by a new 
FSS analysis of MC simulations up to L = 64 along the N line. We obtain p* = 0.89083(3) 
and T* = 0.9527(1), yi = 0.66(1) and y2 = 0.250(2) for the RG dimensions of the two 



^ Indeed, assuming a nondegenerate ground state (this should be the case in two dimensions), we have 
fj,4 = (see App. B for the definitions), from which the relation follows. 



0.09(1). 
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L 


P 


/3min 


Praax 


Nt iVs/64 ATr^n/lO^ 


A^therm/lO^ 


32 


0.8910 


0.800 


1.55 


5 


15625 


800 


240 


32 


0.8915 


0.800 


1.55 


5 


15625 


800 


240 


32 


0.8920 


0.800 


1.55 


5 


15625 


800 


240 


32 


0.8925 


0.800 


1.55 


5 


15625 


800 


240 


32 


0.8930 


0.800 


1.55 


5 


15625 


800 


400 


32 


0.8915 


0.800 


2 


6 


15625 


800 


240 


32 


0.8920 


0.800 


2 


6 


15625 


800 


240 


32 


0.8925 


0.800 


2 


6 


15625 


800 


240 


32 


0.8930 


0.800 


2 


6 


15625 


800 


240 


48 


0.8915 


0.740 


2 


9 


31250 


2000 


400 


48 


0.8920 


0.740 


2 


9 


31250 


2000 


400 


48 


0.8925 


0.740 


2 


9 


31250 


2000 


400 


48 


0.8930 


0.740 


2 


9 


31250 


2000 


600 


64 


0.8915 


0.710 


2 


13 


7813 


3000 


900 


64 


0.8920 


0.710 


2 


13 


7813 


3000 


900 


64 


0.8925 


0.710 


2 


13 


7813 


3000 


900 


64 


0.8930 


0.710 


2 


13 


7813 


3000 


900 



TABLE VI: Parameters of the random-exchange MC runs for L > 32. 

relevant operators in the absence of external field, and rj — 0.177(2) for the magnetic critical 
exponent associated with the spin-spin correlation function. 

In Fig. 4 we report the available estimates of the critical points and report simple in- 
terpolations, discussed in App. C, which take into account all theoretical predictions and 
numerical results. 

Discussions with Marco Picco and correspondence with Nihat Berker and Alexander 
Hartmann are gratefully acknowledged. 
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APPENDIX A: SOME DETAILS ON THE MONTE CARLO SIMULATIONS 



In our parallel-tempering simulations we consider Nt systems at the same value of p and 
at Nt different inverse temperatures /Jmin = Pi, . . . , (3^^ = /Smax, where /Smax is chosen to be 
either 2 or 1.55. To avoid repeating the runs twice, for L — AS and 64, P^aax is always chosen 
to be 2, while one of the Pi corresponds to 1.55. Moreover, for all values of L, we choose 
Pi — Pn{p) for some i, where Pn{p) is given in Eq. (4), so that the corresponding point lies 
on the N line. This choice gives us estimates along the N line, which can be compared with 
exact and previous numerical results. They provide a check of the numerical simulations 
and allow us to improve the estimates of the critical parameters of Ref. [15], see App. C. 

The elementary unit of the algorithm consists in A^ex = 20 Metropohs sweeps for each 
configuration followed by an exchange move. We consider all pairs of configurations cor- 
responding to nearby temperatures and propose a temperature exchange with acceptance 
probability 

V = exp{(A - Pi+i){Ei - Ei+i)}, (Al) 

where is the energy of the system at inverse temperature We generate A^^ disorder 
samples, and for every sample we perform a MC run of Aj-un Metropolis sweeps for each 
Pi value. The first Atherm iterations are discarded for thermalization (see Ref. [15] for a 
discussion of the thermalization issues). The parameters of the runs with L > 32 are 
reported in Table VI. Finally, note that the determination of U22 requires the computation 
of a disorder average of the square of a thermal average. We use an essentially bias-free 
estimator discussed in Ref. [54]. 



APPENDIX B: DEFINITIONS 



The two-point correlation function is defined as 

G{x) = [{a^a,% (Bl) 

where the angular and the square brackets indicate the thermal average and the quenched 
average over disorder, respectively. We define the magnetic susceptibility x = Ylix ^(^) 
the correlation length ^, 



^min^(?min) 
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where gmin = (27r/L, 0), q = 2 sing/2, and G{q) is the Fourier transform of G{x). We also 
consider the magnetization m defined as 

where V is the volume, and the specific heat Cy 

where is the Hamiltonian. 

We also consider quantities (we call them renormalized couplings) that are invariant under 
RG transformations in the critical limit. Beside the ratio 

= i/L, (B3) 

we consider the RG invariant quantities 

J J _ [/^4] _ [A^j] - [^^2f _ J 

where 

Hk^iiJ^^.)'). (B4) 

X 

APPENDIX C: CRITICAL EXPONENTS AT THE MULTICRITICAL POINT 

In each parallel-tempering simulation we fixed p and considered several values of (3 from 
Pmin < P* up to /^max which is either 2 or 1.55, hence larger than the multicritical value 
(3*. In all runs we were careful to include a point on the N line. Since the energy is known 
exactly on this line, this choice allowed us to test the correctness of the simulation code. 
Moreover, we were able to collect a significant amount of new data, which can be combined 
with the old ones presented in Ref. [15]. As we shall see, the FSS analyses of this new set 
of data allows us to improve the estimates of the critical parameters. 

As in Ref. [15] we perform combined fits of the renormalized couplings to Eq. (11) and 
(12). The new results are reported in Tables VII and VIII. The estimates of yi are quite 
stable and essentially independent of L^in, of the observable, and of the scaling corrections. 
We thus quote 

yi = 0.66(1), (CI) 
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-^min 


xVdof 


P* 


yi 


R^,U4,U22 12 


383/289 


0.890864(4) 


0.659(2) 


16 


207/220 


0.890844(4) 


0.658(2) 


24 


120/151 


0.890828(6) 


0.658(3) 


32 


58/82 


0.890822(8) 


0.651(5) 


R^,U^,Ud 12 


424/289 


0.890853(3) 


0.660(1) 


16 


248/220 


0.890856(4) 


0.660(2) 


24 


194/151 


0.890850(5) 


0.659(5) 


32 


100/82 


0.890848(7) 


0.653(5) 



TABLE VII: Estimates of p* and yi at the MNP. Results from combined fits of three different 
renormalized couplings to Eq. (11) with Timax = 2. Here Ua = U4 - U22. 

where the error is chosen quite conservatively, and is such to include all results. This result 
is fully consistent with the estimate yi = 0.655(15) of Ref. [15]. The estimates of p* vary 
between 0.89081 and 0.89086, so that we quote 

p* = 0.89083(3). (C2) 

This estimate agrees with that we obtained in Ref. [15], i.e. p* = 0.89081(7). Moreover, it 
is in full agreement with the recent calculations of Ref. [14] : Two different approximations 
gave p* fti 0.890822 and p* ^ 0.890813. 

Our analyses also provide estimates of the critical-point value of the renormalized cou- 
plings: 

Rl = 0.997(1), (C3) 
[/; = 1.1264(4), (C4) 
U;^ = 0.0817(3). (C5) 

Scaling corrections are particularly weak and apparently decay as or faster. Note that 
this does not necessarily imply the presence of nonanalytic corrections associated with RG 
irrelevant operators with lu ^ 2. Indeed, in all cases we expect contributions due to the 
regular part of the free energy, which decay as L''"^ L~^^^. 

The critical exponent y2 is derived from the critical behavior of R' = dR/dp, where R is 
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i^min X7D0F P* yi U 

R^,U4,U22 6 283/318 0.890822(7) 0.665(3) 1.79(13) 

8 232/300 0.890814(9) 0.660(10) 1.98(25) 

R^,U4,Ud 6 396/318 0.890864(3) 0.662(2) 3.18(10) 

8 332/300 0.890857(1) 0.660(2) 4.11(27) 

TABLE VIII: Estimates of p*, yi, and lo at the MNP. Results from combined fits of three different 
renormalized couphngs to Eq. (12) with nmax = 2 and femax = 1- 

a renormalized coupling [15]. Neglecting scaling correction, its FSS behavior is given by 

R' = ^Ly^h{uiLy\u2Ly^) + ^Ly^h{u,Ly\u2Ly^), (ce) 

where u\ and U2 are the nonlinear scaling fields associated with the two leading relevant 
operators. In general, we expect [36] U2 to vanish on the N line, so that 

U2{(3,p)^S{(3-(3m{p).P-P*). (C7) 

where I3n{p) — ^ITn{p), Tn{p) is defined in Eq. (4), and the function S{x,y) is such that 
S{0,y) — and dS{0,0)/dx ^ 0. Since the transition fines must be tangent to the fine 
p = p* as a consequence of a general rigorous inequality [13], we also have 

p) — p — p* + quadratic terms. (C8) 

The independence of ui on /3 at leading order, implies that the first term in Eq. (C6) vanishes 
at the MNP, so that R' ~ Ly^ for L ^ oo at p = p*. 

In order to compute 1/2, we perform three different fits of our data on the N line. In the 

first one, we neglect the p dependence of du2/dj3 and set dui/djS = 0. Then, setting U2 = 
and expanding in powers of tiiL^^ ^ {p — p*)Ly^, we obtain 

^max 

InR' = y2lnL + Y^ Gnip - p*YU'y^ . (C9) 

n=0 

In the second fit we include the nontrivial dependence of U2 on j3 and p. We fit the results 
to 

In i?' = y2 In L + J] a„(p - p*)"L"^i + ^ &^(p _ p*)-. (CIO) 



n=0 m=l 



24 



-^min 


xVdof 


2/2 




y2 


R'^ 8 


90/102 


0.2533(6) [5] 


57/101 


0.2521(5) [5] 


12 


88/96 


0.2535(7)[5] 


56/95 


0.2519(6) [6] 


16 


60/73 


0.2530(9) [6] 


45/72 


0.2514(8) [7] 


24 


37/50 


0.2531 (13) [7] 


31/49 


0.2515(13) [8] 


32 


16/27 


0.2528(21)[8] 


15/26 


0.2523(21) [8] 


Ui 8 


158/102 0.2492(9) [18] 


90/101 


0.2480(4) [17] 


12 


148/96 


0.2496(11)[19] 


89/95 


0.2478(5) [18] 


16 


95/73 


0.2450(14) [21] 


65/72 


0.2480(7) [21] 


24 


53/50 


0.2509(18) [25] 


42/49 


0.2490(11) [24] 


32 


32/27 


0.2514(23) [28] 


29/26 


0.2500(19) [27] 



TABLE IX: Estimates of 7/2 at the MNP. We fix yi = 0.66(1) and Pc = 0.89083(3). Results from 
fits of the derivative of the renormalized couplings R^ and U4. On the left the results refer to the 
fit to Eq. (C9) with rimax = 2, on the right to the fit to Eq. (CIO) with nmax = 2 and mmax = 1- 
The error in parentheses is the sum of the statistical error and of the error due to uncertainty of 
yi; the error in brackets gives the variation of the estimate as Pc varies by one error bar. 

Finally, note that Ui may depend on /? at quadratic and higher orders, so that on the N line 
one may have 

rill, 

p - p* + 0[{p - p*y]. (Cll) 



dp 

Hence, the first term in Eq. (C6) may give rise to corrections of order [p — p*)V'^~'^'^ . Thus, 
we also perform fits to 

Ini?' = yalnL + ^ a„(p - p*)"L"^i + L'^^ ^ hk{p - p^fh^^K (C12) 

n=0 fe=l 

In Table IX we report the results of the fits of R'^ and U'^ to Eqs. (C9) and (CIO). The 
inclusion of the analjd;ic corrections significantly reduces the and changes sHghtly the 

estimates of y2- Fits to Eq. (C12) give results which are essentially equivalent to those 
obtained by fitting to Eq. (CIO). Comparing all results we obtain the estimate 

y2 = 0.250(2), (C13) 

which is identical to that reported in Ref. [15]. 
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InZ 


Inx 


-^min 


xVdof 






8 


393/102 0.1736(6)[12] 


2653/102 0.1752(2)[5] 


12 


272/96 


0.1747(7) [13] 


2340/96 0.1749(3) [5] 


16 


146/73 


0.1760(8) [14] 


1342/73 0.1751(3)[5] 


24 


68/50 


0.1776(10)[16] 


761/50 0.1752(4)[6] 


32 


30/27 


0.1782(12) [18] 


110/27 0.1761(4)[7] 



TABLE X: Estimates of r? at the MNP. We fix yi = 0.66(1) and Pc = 0.89083(3). Results from 
fits of In Z and In % to Eq. (22) with rimax = 2 without analytic correction (in the case of In Z the 
coefficient of InL is of course —rj). The error in parentheses is the sum of the statistical error and 
of the error due to uncertainty of yi ; the error in brackets gives the variation of the estimate as Pc 
varies by one error bar. 

Finally, we determine rj. We compute it from the critical behavior of x and, as in Ref. [15], 
from that of Z = x/'C^- The results of the fits with and without analytic corrections are 
reported in Tables X and XI. The most stable results are obtained from fits of x which take 
into account the analytic corrections. As final result we quote 

77 = 0.177(2), (C14) 

where the error is such to include the estimates of 1] obtained from the analysis of InZ. This 
result is consistent with the estimate r] = 0.180(5) reported in Ref. [15], but significantly 
more precise. 

The results obtained here allow us to predict the behavior of the different transition lines 
close to the MNP. Standard scaling arguments predict that, close to the MNP, the transition 
lines are given by 

ui\u2\-'^ = X±, (C15) 

where and X_ are two constants that refer to the lines which satisfy T > T* and T <T*, 
respectively. They can be determined by considering the estimates of the critical points Pc, Tc 
close to the MNP. The crossover exponent is equal to the ratio yi/y2- In the present case 
we have 

0=^ = 2.64(5). (C16) 
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InZ 






-^min 




V 


xVdof 


V 


8 


300/101 0.1745(3)[11] 


83/101 


0.1767(1)[5] 


12 


117/95 


0.1763(3) [12] 


79/95 


0.1768(1) [6] 


16 


60/72 


0.1776(4) [14] 


50/72 


0.1774(1) [6] 


24 


32/49 


0.1791(5)[15] 


36/49 


0.1771(1)[7] 


32 


15/26 


0.1794(9) [17] 


20/26 


0.1771(2)[7] 



TABLE XI: Estimates of ?? at the MNP. We fix yi = 0.66(1) and Pc = 0.89083(3). Results from 
fits of In Z and In x to Eq. (22) with 

'^max — 2 and TTimax — 1- The error in parentheses is the sum 
of the statistical error and of the error due to uncertainty of yi; the error in brackets gives the 
variation of the estimate as Pc varies by one error bar. 

We can use Eq. (C15) to obtain an interpolation of our results up to T = 0, which represents 
our best guess of the transition line, given the estimates of the critical points we have. For 
this purpose, we choose 

U2{p, T) = tanh(l/T) - 2p + 1, (C17) 

so that M2 = along the N hne, cf. Eq. (4). Thus, the critical line is given by the approximate 
expression 

Pc - P* + a2{T, - T*f = X_U2{pc, Tcf, (C18) 

where we have kept the 0{AT'^) in the analytic quadratic corrections to the linear behavior 
of Ui, cf. Eq. (C8). Since 2 < < 3 this quadratic term is dominant in the asymptotic 
expansion at the MNP, while the nonanalytic term in the right-hand side of Eq. (C18) 
represents a next-to-leading contribution. Since Pc — P* ~ {Tc — T*y the other quadratic 
terms appearing in the expansion of Ui are subleading. The free parameters 02 and X_ 
are fixed by requiring the hne to go through the points {pc — 0.8925(1), Tc = 0.5) and 
(po = 0.897, T = 0). We obtain 02 = -0.0061 and X_ = 0.0386. The corresponding hne is 
reported (dashed hne) in Fig. 4. The interpolation (C18) gives Pc = 0.89159 at P = 1.55, 
and the derivative dpc/dj3 = 0.00180 ai (3 = 2 which are in good agreement with the MC 
estimates Pc = 0.8915(2) at /5 = 1.55, and dpJdjS = 0.0020(3) at /3 = 2 obtained in Sec. III. 
We have also determined an interpolation of the available numerical data [11, 22] valid 
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for T > T*. A simple expression, which satisfies Eqs. (C15) and (3), is 

p^ = p* + {13* - /3)^ '^^(1.41484 - 4.25764/5 + 5.67965/3^ - 2.77095/?^), (C19) 

with p* — 0.89083 and P* — 1.04962. The corresponding fine is reported (dashed fine) in 
Fig. 4. 
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